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Abstract 

The Bogoliubov equations are solved for a three-dimensional Bose-Einstein condensate containing 
a doubly quantized vortex, trapped in a harmonic potential. Complex frequencies, signifying 
dynamical instability, are found for certain ranges of parameter values. The existence of alternating 
windows of stability and instability, respectively, is explained qualitatively and quantitatively using 
variational calculus and direct numerical solution. It is seen that the windows of stability are much 
smaller for a cigar shaped condensate than for a pancake shaped one, which is consistent with the 
findings of recent experiments. 
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I. INTRODUCTION 



One of the hallmarks of superffuid flow is the quantization of fluid circulation. Its mani- 
festation through the formation of quantized vortices is a long studied subject in supercon- 
ductors, superfluid Helium and trapped atomic gases Although a vortex can in principle 
carry any number of circulation quanta, it is well known that in an infinite homogeneous 
system, only singly quantized vortices are energetically allowed. In such a geometry, a vortex 
with quantum number larger than unity, frequently called a giant vortex, has higher energy 
than the corresponding number of separated vortices with single quanta while carrying the 
same angular momentum. 

For Bose-Einstein condensates trapped in magnetic potentials, the situation is more com- 
plicated since the size of a vortex core can be comparable to the total system size. Still, it 
is an accidental fact that quantum numbers larger than unity are energetically unfavorable 
in Bose-Einstein condensates in parabolic potentials, just as in the homogeneous case. This 
result hinges on the fact that in the noninteracting limit, the lowest energy levels with a 
given angular momentum are vastly degenerate; interactions lift the degeneracy and favor 
a state with a diluted density profile, which is a lattice of singly quantized vortices. The 
situation is different if the trapping potential is steeper than harmonic, in which case giant 
vortices are energetically favorable in the limit of weak interactions 

The related but separate issue of dynamical stability of doubly quantized vortices has 
become especially relevant lately since the problem can be readily probed by experimental 
means jjj. Theoretically, it is known that in the two-dimensional limit such a vortex is 
dynamically unstable towards splitting into two only in certain windows of parameter space, 
when the interaction strength lies within certain intervals [4]. These intervals have been 
shown to be bounded by zeros of eigenvalues for the corresponding static stability problem 
jjj. However, the situation is entirely different in three dimensions: in the strongly cigar 
shaped limit, instabilities are expected to lie densely in phase space, because of the closely 



spaced energy levels in the third direction jo, 13, |S( . In this paper, we perform a systematic 
study of the Bogoliubov excitations in two- and three-dimensional condensates containing a 
doubly quantized vortex. The paper is organized as follows. In Sec.|H]we discuss the equa- 
tions and the structure of the stability problem. In Sec. IHII we address the two-dimensional 
problem and show how the numerically observed parameter dependence can be understood 
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and approximated analytically. Section IIVI contains an overview of the numerical results 
for the three-dimensional condensate. In Sec. IS we discuss in greater detail the limit of a 
cigar shaped condensate. Section |VU relates the present results to experimental findings. 
We calculate the instability regions for an anharmonic trap in section IVm to investigate the 
stability in a more general case. Finally, in Sec. IVllJ we summarize and conclude. 

II. STABILITY PROBLEM FOR A CONDENSED BOSE GAS 

The system under study is a Bose-condensed atomic gas that is trapped in a cylindrically 
symmetric harmonic potential. At zero temperature in the dilute limit it is described by a 
condensate wavefunction \l/(r, 9, z, t) that obeys the Gross-Pitaevskii equation 



The eigenvalue [i is the chemical potential, and the wavefunction \1/ is normalized to the 
number of atoms N. The trapping potential is assumed to be of the form 



The anisotropy of the trap is governed by the ratio A between the axial and radial trapping 
strengths, so that when A is larger than unity, the atomic cloud resembles a pancake and 
when A < 1, it is cigar shaped. We shall see that the system behaves very differently in 
these two limits. The inter-particle interactions are parametrized by an s-wave scattering 
length a, so that Uq = 47ih 2 a/m. Combining these parameters into dimensionless quantities, 
we see that the physics of the system is entirely determined by A and the effective coupling 
strength C = 4nNa/a osc , where a osc is the harmonic oscillator length, a 2 osc = h/(mu)±). 
Since the chemical potential /z is a monotonically increasing function of C (although explicit 
expressions can be found only in limiting cases), it is possible to switch to the pair of 
parameters (A, /i) instead of (A, C). This parameter choice is found to be more helpful in 
explaining the physics of the problem and is therefore used in most of this paper. We now 
switch to dimensionless units in which energy is measured in units of fiiu^, frequency in units 
of u± and length in units of a osc . The Gross-Pitaevskii equation becomes 



r n 2 



V 2 + V(r) + U \^\ 2 * = /i*. 



(1) 



2m 




(2) 




(3) 
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We shall, however, find it natural to reinsert units when discussing physical scales throughout 
the paper. 

For a given stationary solution of the Gross-Pitaevskii equation with eigenvalue /x, 
the small-amplitude excitations of the condensate are defined through the ansatz 



tt(r,t) 



e - *"*, (4) 



*o(r) + K(r)e-^* + ^(r)*e^) 

n 

where u n and u n are the eigenvectors and u) n the eigenvalues of the Bogoliubov equations, 

B(p n (r) = u n (p n (r), (5) 
where <f n { r ) = (u n ( r ), v n( r )) T i an d the linear operator B is defined by 

B=( + V(r) ~ fI + 2C ^ 2 C * 2 \ (6) 

\ -C(^*) 2 -(-|V 2 + V r (r)-/i + 2C|^| 2 ) J ' 

B is non-Hermitian and may have complex eigenvalues. If this is the case, the system is 
dynamically unstable and the corresponding modes will grow exponentially, as seen from 
Eq. (j3J). For the problem at hand, where \I/ is assumed to contain a doubly quantized 
vortex at the origin, it is known in the two-dimensional limit that there exist intervals of the 
coupling constant C where exactly one pair of eigenvalues is complex; in these intervals the 

doubly quantized vortex is unstable towards splitting, and between them it is dynamically 

rn 

stable. On the other hand, in the strongly cigar shaped limit, the experiment of Ref. p| saw 
instability for all coupling strengths. In this paper, we shall see how these two situations 
emerge as the extreme limits of the general three-dimensional problem. 

We write the matrix element of a two-by-two, possibly space-dependent, matrix A be- 
tween two Bogoliubov eigenvectors as (<p m \A\<p n ) = J d 3 r(p m (ry A(r){p n (r) . The "inner 
product" of two Bogoliubov eigenvectors ip n is defined with the help of the Pauli matrix 
o z = diag(l, —1) as 

(<Pm\crz\<Pn) = J d 3 r (w m (r)*M n (r) - v m (r)*v n {r)) . (7) 

The "norm" of a Bogoliubov eigenvector <p n with respect to this "inner product", (fnWzl^Pn) , 
can be either positive, negative or zero, since it obeys the relation 

{uj n - u* n )(ip n \cr z \ip n ) = 0. (8) 
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If u n is complex, the norm must be zero. In the case of real u n , the norm may be either 
positive or negative and we shall fix normalization such that its absolute value is unity for all 
modes with real eigenvalue. Note that for every solution ip n = (w n ,t>*) T of the Bogoliubov 
equation with eigenenergy u n , there exists a solution (p n = (t>*,-u n ) T with eigenvalue —a;*. 
If \l/ is the ground state of the Gross-Pitaevskii equation in the absence of rotation, the 
norm of each mode <p n is equal to the sign of the corresponding u n , but that result does not 
hold for a general The existence of modes with positive norm and negative energ 



y 


'and 


9, 


12 1 . 



If the problem is cylindrically symmetric and the condensate has an angular momentum 
M per particle, the excitations can be labeled by an angular momentum quantum number 
m relative to that of the condensate, such that 



u{r,9,z) = e l ^ M+m)0 u(r, z), 



J,(M+m)$ 

v{r,9,z) = e iiM ~ m)9 v(r, z). 



(9) 



The Bogoliubov eigenvalue problem thus splits up into a block diagonal matrix where the 
blocks corresponding to different m are decoupled. This allows us to treat each m value 
separately; the m'th sector of the Bogoliubov matrix takes on the form 



B 



D x + V d {r,z) V c {r,z) 
-V c (r,z) -(D 2 + V d (r,z)) 



(10) 



where 



D 2 



d 2 1 d 
dr 2 r dr 
d 2 1 d 



(M + m) 2 | d 2 



(M - m) 



+ 



dz 2 
dz 2 



(11) 



and the diagonal and off-diagonal potentials are defined as 



V d (r,z) = V(r,z) - fx + 2C\^(r,z) 



V c (r,z) = CV(r,z) 2 . 



(12) 



For the problem at hand, the m = 2 sector will turn out to be especially interesting, as we 
shall see in the following sections. 

Let us now study what happens when two Bogoliubov eigenvalues collide as a control 
parameter is changed. We will see that they either become complex, signifying a dynamical 



instability, or undergo an avoided crossing. Consider two Bogoliubov amplitudes ipj, with 
j = 1,2, that solve the Bogoliubov equations for some pair of parameter values A, C. We 
assume that the corresponding eigenvalues uij are real. Now slightly increase C to a value 
C = C + 5C, so that the ground-state solution \I/ becomes = ty + S^fr, and the Bogoliubov 
matrix is correspondingly written B' = B + 5B. Consider the restricted Bogoliubov problem 
in the truncated basis formed by ifix,f2- (The validity of this approach shall shortly be 
determined.) Writing the ansatz as cp — ax(fx + a 2 ip 2 , one obtains the eigenvalue problem 
for the vector a = (ax, a 2 ) T , 

(<PiWz\<Pi) \ I Ux + SBxx SB 12 

(f2Wz\f2) I \SB 21 UJ 2 + 5B 22 

where SB^ = (ipi\a z SB\ipj). Since a z SB is a Hermitian matrix, we have SB 2 x = 5B 12 * . The 
eigenvalues are 

ujx + 5B X x + oj 2 + 5B 22 
U± = 2 

2 



a = una. 



(13) 



± 



1/2 

ui-u 2 + 6B n -6B&Y |2 

1 ' Wl\(7zWl)W2\(7zW2)\0B 12 \ 



(14) 



2 

It is seen that the eigenvalues can be substantially altered only if the energy difference ux —u 2 
is of the same order of magnitude as the perturbations SBij; thus one only has to care about 
the coupling between nearby energy levels, so the truncation of the basis to two states 
is justified. When the two modes have equal norm, they experience an avoided crossing. 
Conversely, a condition for the eigenvalues to be complex is that the two modes have different 
norm, (v?i|^|v 9 i)( ( / 3 2|^|¥ , 2) = — 1- Thus it can be concluded that a dynamical instability 
is formed by the mixing of a Bogoliubov mode of positive norm with one of ne gati ve norm 
(states with opposite Krein signature in the language of Hamiltonian systems [16]); when 
their energies coincide they go over to a pair of zero-norm modes with eigenvalues that are 
complex conjugates of each other. We will see plenty of examples of this phenomenon in the 
following. 

For the numerical solution of the Bogoliubov equations, we employ the following method. 
At each point in the (A, /i) phase space, the stationary solution of the Gross-Pitaevskii 
equation ((TJ) is sought under the condition that it contains a doubly quantized vortex along 
the z axis, so that the wave function is assumed to be of the form 

V(r,6,z) = e 2ie ty(r,z). (15) 



After that, the Bogoliubov problem is a matter of finding the eigenvalues to a 2n x 2n 



matrix, where n is the size of the numerica 



discrete- variable representation (DVR) 



aid. The spatial coordinates are discretized in a 
using a Laguerre mesh for the radial direction 



and a Hermite mesh for the axial direction. For the solution of the two-dimensional problem 
reported in Sees. 11111 and IV111 the z dependence drops out and a one-dimensional Laguerre 
mesh is employed. 

III. TWO-DIMENSIONAL CASE REVISITED 

We first consider the case of a two-dimensional condensate, which can be realized in a 
strongly pancake-shaped trap, i. e. in the limit of large A. In this limit, the condensate is in 
the harmonic-oscillator ground state in the z direction; integrating out the z-dependence, 
one obtains the effective coupling constant C = CX 1 ^ 2 / V2~rr. As mentioned above, we choose 
to state our results in terms of the chemical potential /i, which is a monotone function of C; 
this will simplify some of the calculations and make some relations more clear. 

In the two-dimensional case it is known that the doubly quantized vortex is dynamically 
unstable only in certain intervals of the coupling strength [4] . The first instability sets in at 
zero coupling, as can be seen by considering the restricted eigenvalue problem, Eq. (fT^j) . for 
the lowest-lying eigenstates of the noninteracting problem. Thus we take <pi = (04oo, 0) T and 
<£>2 = (0, 0ooo) T , where (pi, nr ,n z (r,9, z) are the harmonic-oscillator eigenstates in cylindrical 
coordinates, 



</W 2 (r,M = J^l— H nz (z)e~ Xz2/ \ T^L^r 2 ) {re^e^l\ (16) 
where H n is a Hermite polynomial and 

^ \n-jjj\ 

is a generalized Laguerre polynomial. The two eigenstates <fi and ip2 are degenerate with 
Bogoliubov eigenvalue uj = 2 in the noninteracting limit, but have opposite norms. Thus, 
for any finite coupling strength C, the Bogoliubov eigenvalues calculated from Eq. (fTHjl turn 
complex. In order to determine the point where the system turns dynamically stable again, 
one has to consider the contributions from states with higher n r j^|. 



7 



Figure ^ shows the numerically computed real parts of the Bogoliubov eigenenergies 
as functions of the chemical potential \i in two dimensions. According to Eq. (jl4j) . the 
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FIG. 1: Energy levels in two dimensions for m = 2. The * represent the analytical approximation 
to the crossing of the core mode with the trap modes shown in Eq. ()2U|) . 

merging of two lines into one signifies the onset of instability, as two real energies become 
a complex conjugate pair. It is seen how the successive instability windows arise from 
the coupling of one mode, whose energy increases monotonously, with successive radially 
excited modes. This mode has no radial nodes and has negative norm in the intervals 
where its energy is purely real; we refer to it as the core mode since it is confined to the 
core of the doubly quantized vortex and causes a splitting of the vortex into two. Its 
physical significance is that the doubly quantized vortex is energetically unstable towards 
splitting, which furthermore causes dynamical instability in certain windows of parameter 
space. As the coupling increases, the energy of the vortex is negligible compared with 
the energy contributions from the interaction and trapping potential. Therefore all the 
Bogoliubov energy levels approach the excitation energies for a two-dimensional condensate 
in the Thomas- Fermi limit ^| , except for the core mode whose (real part of its) energy 
increases with increasing fi. In fact, the dependence of the core mode energy on the chemical 
potential is approximately linear. This can be understood intuitively as follows. If the core 
mode is populated, the doubly quantized vortex is split into two. The energy of the core 
mode corresponds to the frequency of the precession of the two vortices (since it determines 
the rate of change of the relative phase of the condensate wave function and the Bogoliubov 
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mode). Each of the two vortices moves with the local velocity, which is given by the velocity 
field from the other vortex and is proportional to the inverse of their separation [lj]. The 
angular frequency of this precession is thus d; pre c = V( mr i2)- Inserting for r 12 the healing 
length £, which is the size of a vortex core, one obtains indeed hu) prec oc [i 

A more quantitative estimate of the core mode energy can be obtained in the limit of 
a large condensate. Consider the Bogoliubov equation (jSJ) with Eq. (fTUj) inserted for the 
matrix operator B, and neglect the z depencence. We do the analysis for a general set 
of quantum numbers M, m before specializing to the currently relevant case M = m = 2. 
The diagonal and off-diagonal effective potentials Vd(r) and V c (r) defined in Eq. (fTTH) are 
sketched in Fig. |21 Since the norm of the core mode ip is negative, it is dominated by its 
lower index v. The lowest-lying eigenstates of the diagonal potential Vd are concentrated 
to the vortex core, i. e. the innermost potential well of the potential Vd- This justifies 



FIG. 2: Effective potentials entering the Bogoliubov equation for a doubly quantized vortex. Full 
line represents the diagonal potential Vd, the dashed line represents the diagonal potential with 
the addition of a centrifugal term, Vd + 4/r 2 , and the dotted line is the off-diagonal potential V c . 
The potentials are defined in Eq. (|10|). 

the term "core mode" for this mode, and we therefore denote the amplitude by t> core and 
the corresponding energy by u; C orc- The corresponding upper amplitude w core experiences 
an additional centrifugal potential which pushes it away from the origin. As a result, the 
overlap between u core and t> core weighted by V c is exceedingly small and we can neglect the 
off-diagonal term in the Bogoliubov equation for this mode. Moreover, in the limit of a 
large condensate, the extent of the function f CO re( r ) is of order £ and the trapping potential 
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contributes little to its eigenenergy. Hence the Bogoliubov equation for the core mode can 
be well approximated by the diagonal term neglecting the external potential, 



h 2 ( d 2 Id (M — m> 2 



-7T + " 5-^- + 2 W(r 



-uJ COTe v core {r). (18) 



2m \ dr 2 r dr r 2 
The condensate wave function \I/ can be written as the product of the square root of the 
density away from the core, n , and a core function f(r) describing how \l/ goes to zero due 
to the centrifugal force: \l/(r) = y/riof(r). The density no is unaffected by the trapping 
potential at distances r comparable to £, and is assumed constant; it is related to the 
chemical potential by fi = U n . Rescaling the length, r = £x, results in the equation 



1 / d 2 Id (M-m) 2 \ 1 p . 



2\ dx 2 xdx ' x 2 J 2 ^ v-/ j -orc(x) - ™ c v core (x). (19) 
Already from here one can see that since all parameters are scaled away from the left-hand 
side, u core /fi must be a constant, i. e. the core mode energy is proportional to the chemical 
potential. A variational quantitative estimate can be obtained as follows. Assuming a 
general condensate with angular momentum M and modeling f(x) as a linearly increasing 



function cut off at the position x = 6, we obtain by variational means b = My/E [ll|. With 
this assumption for the function /, the solutions f core to Eq. (fT9"j) are harmonic-oscillator 
eigenfunctions and the eigenenergies come out as 

^ = -<^ + <"-»•) + 1)^ + 1. PO) 

We are interested in the core modes with n r = 0. This energy is indeed positive for M > 2 
and m > 0, i. e. of opposite sign compared to the norm. (Recall that we could equally 
well have considered the corresponding Bogoliubov mode with m — > — m and the upper 
component u confined to the core; it has positive norm and negative energy.) 

The above estimate is readily checked against the numerical result for a two-dimensional 
trapped condensate in the strong-coupling regime. We read off the numerically calculated 
core mode energy at the arbitrarily chosen point \x = 40fkj±. For the m = 2 core mode 
in a condensate with M = 2, we obtain the variational estimate u> cor e = 0.423/x, while 
the numerical result is a; core = 0.438/x. For m = 3 and M = 3, we obtain variationally 
^core = 0.615/i and numerically u> core = 0.665//, and for the m = 4 core mode in a M = 4 
condensate the variational result is owe = 0.711// and the numerical result is u cove = 0.781//. 

The positions of the successive instability windows can now be estimated by calculat- 
ing the crossing of the core mode with quadrupole modes with increasing radial quantum 
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numbers n r . The energies of the latter will in the limit of strong coupling not depend on 
the multiply quantized vortex in the center of the condensate, so we can use the values 
calculated for a nonrotating two-dimensional condensate jl^, [l4( 

u nr ,m = ?&j± \j % n r + 272 r m + 2n r + m. (21) 

The variational estimates for the crossings, i. e. the points fi where cu core = u) nr>m for m = 2 
and some n r , are indicated in Fig. [T] It is possible to proceed and calculate the overlap 
between the core and quadrupole modes, which yields the widths of the unstable windows 
and the imaginary parts of the eigenfrequencies. However, this is impossible in practice, 
since the result will be extremely sensitively dependent on the width of the core mode, 
which is a variational parameter. We will therefore not pursue this analysis. 



IV. INSTABILITY REGIONS FOR A THREE-DIMENSIONAL CONDENSATE 

Figure H3 contains the main result of this study. Displayed is the largest imaginary part 
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FIG. 3: Unstable regions for a doubly quantized vortex in the parameter space of trap anisotropy 
A and shifted chemical potential fj, s . The shading indicates the largest of the imaginary parts of 
the Bogoliubov eigenvalues. White is zero which means that the condensate is stable. 

of the Bogoliubov eigenenergies as a function of trap anisotropy A and shifted chemical 
potential /x s . The latter is defined as the chemical potential downshifted by its value in the 
noninteracting limit, 

fi a = H-X/2-3, (22) 
11 



and it will play the role of our coupling parameter in the following. All the unstable modes 
are in the m = 2 sector, there are no unstable modes with quantum number m 7^ 2. An 
instability with m = 2 corresponds to a splitting of the doubly quantized vortex line into 
two single vortex lines, as can be seen by considering the density profile of a superposition of 
m = 2, m = and m = 4 radial eigenfunctions: In the noninteracting limit the condensate 
wavefunction at some time instant t can be written 



1/4 r 



-=(x + iy) 2 + r)—=(x + iyf + rf 
V2 V4! 



where 77 is a complex, time-dependent amplitude. The quantity within brackets can for 
small x, y be written as [x — xq + i(y — yo)][x + xq + i(y + yo)], which describes two singly 
quantized vortices at opposite sides of the z axis, with Rer] = y% — x\ and Imi] = 2x y . Note 
that when the coordinate along the abscissa in Fig. El is chosen to be the shifted chemical 
potential fi s , the instability regions in the large- A limit appear as vertical stripes, as we 
shortly discuss. The instability regions bounded by straight diagonal lines are instabilities 
between negative energy states n r = 0,n z = n n and positive energy states n r = 0,n z = n p 
with n n — n p even. They correspond to instabilities with axial nodes. In the next section 
we shall study the instabilities in greater detail. Of course, one has to keep in mind that 
the quantum numbers n r , n z only make sense in the weak coupling limit. Their meaning is 
lost in the strong-coupling limit, where the r and z dependence is no longer separable, but 
it remains a convenient way to label the states. 

The alternating stability and instability regions known from previous two-dimensional 
studies (4,0,!^ and discussed in Sec. lIIII are clearly seen in the pancake shaped, large A limit 
as vertical stripes. We refer to these as 2D instabilities, since they arise from dynamics in 
the plane. When A approaches unity from above, one sees how the 2D instability regions 
become distorted. In fact, the distortion of the second vertical stripe can be explained as 
an avoided crossing phenomenon. As discussed in Sec. IIII[ it is the crossing of the n r = 1, 
n z = mode with the n r = 0, n z = core mode in the pancake shaped limit that gives rise 
to the second vertical stripe. On the other hand, in the cigar shaped limit the instability of 
the n r = 0, n z = 2 mode with the n r = 0, n z = core mode appears as a diagonal band 
in the lower part of Fig. H2 as we discuss in Sec. El The two modes (n r = l,n z = 0) and 
(n r = 0, n z = 2) have the same symmetry and positive norm, and will mix when their energy 
is similar, giving rise to an avoided crossing around A = 1. This avoided crossing of the real 
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FIG. 4: Real parts of the Bogoliubov eigenenergies as functions of the shifted chemical potential 
fj, s for the trap anisotropy A = 1.2, i. e., an almost isotropic trapping potential. 

energy levels is shown in Fig. H]for A = 1.2. The two modes in question start off at /i s = as 
harmonic-oscillator eigenstates with energies uj nr= \ }nz= Q = 4 and uV-=o,n*=2 = 2 + 2A = 4.4, 
respectively. The lines can be seen to undergo an avoided crossing around fi s = 0.5; the 
exact position of this crossing depends on A. Furthermore, both of these eigenmodes become 
unstable when they cross with the core mode, the negative-norm eigenmode n r = 0, n z = 0. 
The result is an avoided crossing of the instability regions. In the second panel of Fig. El the 
instability region in the high A region corresponds to n T — 1, n z — 0; this goes continously 
over to n r = 0, n z = 2 in the lower left corner, while we find the n r — 1, n z — state again 
as the strong instability between 5 < fi s < 6 in the lower part of the figure. 
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V. INSTABILITY IN THE CIGAR SHAPED LIMIT 



When A < 1, the condensate obtains an elongated, cigar-like shape. In this limit the 
lowest-energy excitations are in the z direction, since their energy separation is \ftw± while 
the radial excitation energy is equal to The problem is thus one- dimensional as long 

as radial excitations can be neglected. The instability regions carrying different quantum 
numbers n z are seen to spread out like a fan in the lower part of Fig. El The analysis 
in the previous paragraph still holds for the m = 2,n r = n z = instability, which sets 
in at zero coupling also in the cigar shaped limit. Because of parity, this mode does not 
mix with the modes that have odd axial quantum numbers n z and we now deal with those 
separately. Consider again the Bogoliubov equation restricted to the space spanned by 
two nearby modes, Eq. ([13)1 . Assuming that perturbation theory holds so that harmonic- 
oscillator eigenfunctions can be used, we insert into Eq. (JT3*|) the trial Bogoliubov amplitudes 
<fi = (</>4oi, 0) T and cpx = (0, 0ooi) T ; where we remind the reader that the harmonic-oscillator 
eigenfunctions <p^n r ,n z are indexed with the azimuthal, radial, and axial quantum numbers 
in turn. The calculation shows that the mode is unstable when 

512tt 3 / 2 C 512tt 3 / 2 
17\/2 + 8^ < A < 17^2 - 8^' 

We rephrase this in terms of the shifted chemical potential, which in the weak-coupling limit 
varies as /x s = 3C*/[3(27r) 3 / 2 ]. There results 

_^<^<^1^. (25 ) 
17 + 4^ A 17-4^6 

Thus, the fan-like structure in the lower part of Fig. El should be bounded from the left 
by the ray A ~ 0.6002/i s . Although the boundary of this instability region appears fairly 
straight, its slope is about half the value predicted by this weak-coupling analysis. This is 
because the instability sets in at a coupling strength where many more harmonic-oscillator 
eigenfunctions are already mixed into the relevant Bogoliubov modes. Going to the limit of 
small A and small fi s does not help, since the above analysis indicates that the instability sets 
in only when the interaction energy, which is proportional to the shifted chemical potential 
fi s , is of the order of the the harmonic-oscillator level spacing \hu>±. 

Figure El shows the real parts of the Bogoliubov eigenenergies as functions of the coupling 
for the fixed anisotropy A = 0.2. It is clearly seen that the radially excited states, in the upper 
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FIG. 5: Real part of the Bogoliubov energy spectrum for the case A = 0.2, i. e., in the limit of a 
cigar shaped condensate. 

part of the figure, are well separated from the lowest axially excited states and thus they do 
not affect the dynamics in the weak-coupling limit. Again, the simplest way to determine the 
quantum numbers of a particular mode is to inspect its energy in the noninteracting limit. 
The two modes that in the noninteracting limit have energies u = (2±X)hw± = (2±0.2)hw± 
are the ones with quantum number n z — 1. The initial linear behavior of the energies can be 
obtained from the harmonic-oscillator eigenfunction analysis in the previous paragraph, but 
the curvature of the lines becomes important, and as a result the levels merge at a higher 
chemical potential than predicted by the weak-coupling analysis. 

As discussed above, the merging of the two lines into one signifies the onset of instability, 
where according to Eq. (fhfj) the two real energies become a complex conjugate pair. At a 
somewhat larger [i s the modes with n z = 2 merge and become complex. The levels split 
apart and become purely real again for a larger value of /i s . As long as the system is 
approximately one-dimensional in the sense that radial excitations are well separated from 
the axial ones, the curves presented in Fig. El are universal and are only dilated by a factor 
A as the anisotropy is changed. As a result, the instability boundaries appear as straight 
lines that give rise to the fan-shaped structure in Fig. El 

In the extremely cigar shaped limit, we can completely ignore the excitations in the radial 
direction. In this limit the problem simplifies enough to allow for analytical calculation of 
the mode frequencies; however, as we shall see, the model is only accurate for very weak 
coupling. The Bogoliubov equation will in this case be an eigenvalue problem of two cou- 
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FIG. 6: Imaginary parts of the Bogoliubov eigenvalues, Im(u>i), in the limit of weak interaction and 
strongly cigar-shaped geometry, calculated using the one-dimensional model defined in Eq. 1)26(1 . 

pled one-dimensional differential equations. Inserting the lowest radial harmonic-oscillator 
eigenfunctions for u and v in the Bogoliubov matrix operator, Eq. (|1U|). integrating out 



the radial direction and redefining units by putting z = (Xz) / \/2/j^ and the axial condensate 
function ^(z) = [i s ty(z)/C we obtain 
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(26) 



If we furthermore assume that the Thomas-Fermi (TF) approximation Q for the condensate 
wavefunction is valid in the z direction, the peak value of \l/ is just equal to 1, and the only 
remaining parameter in the equation is /i s /A. Thus the points at which the different axial 
modes become complex are given by fi s / A = const, i.e, straight lines in figureE! Furthermore, 
in the limit fi s / A — > oo we can neglect the derivative of the function at the center of the 
trap and the eigenvalues u>i/[i s have to approach the eigenvalues of the matrix 
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(27) 



The result is Ui = /i s (— 1/24 ± iy/47 /24) — 0(A). In the actual three-dimensional situation 
we have seen that the imaginary parts of the frequencies rise to a maximum value and 
then decrease to zero; this behavior is not predicted by the ID model and is thus a three- 
dimensional effect. In order to estimate the maximum values of the imaginary parts of 
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the frequencies, we observe that the radial dynamics is expected to begin to matter when 
fx s exceeds unity, so putting fx s = 1 gives an estimate Im(cjj) ~ 0.3. The actual value 
from the three-dimensional calculation turns out to be 0.14 for the maximum of the first 
imaginary frequency, with a slight dependence on A. We conclude from the one-dimensional 
calculation that the periodic regions of stability and the maximum of the imaginary part of 
the eigenvalue is an effect arising from the dynamics in the radial direction. 

It is again instructive to visualize the instability by considering the shape of the wave 
function ip(r, 6, z) with a small admixture of the Bogoliubov amplitudes u and v, as we did in 
Eq. (J2BJ). One sees that the m = 2, n z = modes simply correspond to a straight splitting of 
the doubly quantized vortex into two. The m = 2,n z = 1 mode corresponds to two vortices 
that split at the edges of the condensate, at large \z\, but are joined at z — 0, thus forming 
an X-shaped structure. In experiments one will generally be in a regime where more than 
one mode is unstable. This will result in a intertwining of two vortices. This was studied 
_ lly in Ret, flflB. The siting w aS tound to nucleate in certain interval of „ 
corresponding in the present picture to a high quantum number n z , as is expected for strong 

n 



coupling. It was also proposed in Refs. [a, |£( that the criterion for local splitting can be 
found from a local-density approximation of sorts, by treating the elongated condensate as 
a stack of two-dimensional slices. If the local density integrated over the x-y plane matches 
the instability criterion for the two-dimensional system, an instability can be nucleated at 
that point. This kind of analysis presumably holds in the limit of a large condensate, where 
local-density approximations are expected to hold. 



VI. EXPERIMENTAL LIFETIME OF A DOUBLY QUANTIZED VORTEX 

In the experiment carried out by Shin et al. j2|, a doubly quantized vortex was topolog- 
ically imprinted in a 23 Na condensate and its subsequent decay was tracked by observing 
time-of-flight density profiles. It was argued in Ref. 7] that the initial occupation of the 
dynamically unstable modes, which as we have seen have quadrupole symmetry, is mainly 
because of gravitational sag in the trap during topological imprinting, which produces a 
quadrupolar deformation. In order to compute the lifetime of a doubly quantized vortex 
as observed in the experiments, one has in principle to model the full dynamical process 
including the initial seeding of the unstable modes, their growth and mixing, nonlinear ef- 
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fects that occur once the mode occupation becomes appreciable, migration of the density 
fluctuations along the vortex axis [8| , and finally the expansion of the atom cloud before ob- 
servation of two separate density depressions. However, it is clear that the main contribution 
to the lifetime is given by the rate of exponential growth of the unstable modes. This rate 
is just the maximum imaginary part of the complex eigenvalues (MCE) of the condensate 
at that particular point in phase space. The lifetime depends only logarithmically on the 
initial mode population and is thus insensitive to the seeding process. Whether nonlinear 
processes during the latter stages of the decay can appreciably affect the dependence on 
coupling strength is more of an open question; this remains to be investigated and we shall 
see in this section that the overall parameter dependence seems to be very well described by 
the sole parameter that is the MCE. 

The experiment by Shin et al. [3] was done with aspect ratios A ranging from 1/100 to 
1/20. This is clearly in the cigar-shaped domain. In this region the unstable modes are 
quadrupole modes with different numbers of axial nodes, as we have seen in the preceding 
section. Figure U\ shows the value of the MCE for A = 0.2, as a function of the effective 
two-dimensional coupling strength an z , which is the parameter that was used in Ref. j^j]. 
It is defined as an z = aj ^(r,0)2irrdr, and has to be computed numerically for each data 
point. Because of numerical limitations, we have used a larger value of A than in the 
experiment, but as we have seen, the most important features, such as the position of the 
2D instabilities, do not change when A is decreased. In the region between the two 2D 
instabilities, the most unstable mode acquires successively higher numbers of axial nodes. 
The instability gets weaker with higher density in the region < an z < 12, which is precisely 
the parameter interval that was scanned in the experiment. Clearly, because of a coincidence 
the experiment was performed in the parameter regime lying exactly between the two first 
2D instabilities. 

As we saw in Sec. El finding the MCE is a purely computational task, since a good 
analytical approximation seems to be difficult to construct. Qualitatively we can understand 
that the MCE decreases with increasing coupling because the coupling matrix element is 
smaller when the unstable modes have axial nodes. To compare with the experiment [3], we 
identify the lifetime of a doubly quantized vortex with the reciprocal of the MCE. The result 
is presented in Fig.|H]in the units used in the experiment. The numerical result has the correct 
qualitative behavior, and the quantitative scale is also similar to the experimental result. We 
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FIG. 7: The maximal imaginary part of the complex eigenvalues for a condensate with a doubly 
quantized vortex in the cigar shaped limit with trap anisotropy A = 0.2, as a function of two- 
dimensional coupling strength. The two large peaks are termed 2D instabilities since they do not 
depend on the coordinate along the vortex line, while the unstable modes connected with the 
smaller peaks have nodes in the axial direction. 

conclude that while nonlinear effects and mode mixing may be necessary to quantitatively 
fine-tune the splitting times, it seems that simply taking the reciprocal of the MCE is 
sufficient to obtain the parameter dependence of the vortex lifetime both qualitatively and 
roughly quantitatively. 



VII. STABILITY IN AN ANHARMONIC TRAP 



In the previous sections we have studied a doubly quantized vortex in a harmonic trap. 
In this case the vortex becomes unstable already for any finite value of the coupling C . We 
have seen that this is due to the degeneracy of the eigenvalues in the harmonic trap. To 
understand the connection between the energy spectrum and the instability further, we now 
present the corresponding calculations for an anharmonic trap. We use the potential 

V(r) = -r 2 + (ar f . (28) 

The regions of instability for this potential are shown in Fig. El The energies in the 
noninteracting limit are no longer degenerate as they were in the harmonic case: For a given 
anisotropy, the energy of the negative norm state is below the energy of the positive norm 
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FIG. 8: Time scale for splitting of a doubly quantized vortex in a cigar shaped condensate with 
anisotropy A = 0.2, defined as the reciprocal of the largest imaginary part of the Bogoliubov 
eigenvalues. Units on the axes are chosen to correspond with the experiment of Ref. 

state for weak coupling. The nonlinear coupling therefore needs to attain a finite value in 
order for the two real eigenvalues to meet, so that the eigenvalues turn complex and the 
vortex becomes unstable. After this point the situation will be similar to the harmonic 
case, except that the mode frequencies are shifted, which causes a corresponding shift in the 
positions of the unstable regions. 

In the weak-coupling region, before the first complex eigenvalue appears, it is possible to 
find a rotating frame with rotation frequency Q such that all the positive-norm and negative- 
norm states have positive and negative energy respectively (i. e., all modes have a positive 
Krein signature [16(). If this is the case for all sectors with different angular momentum 
quantum number m (which, in fact, it is for certain Q and C values), then according to 

is energetically stable in 

that rotating frame. This certainly implies that the vortex is stable (both linearly and 
nonlinearly) , both to perturbations of the vortex state and also to perturbations of the 
Hamiltonian that are rotationally symmetric. 
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FIG. 9: Real part of the energy spectrum for the two-dimensional case in an anharmonic potential 
with anharmonicity parameter a = 0.92. 
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FIG. 10: Imaginary part of the energy spectrum for the two-dimensional case in an anharmonic 
potential with a = 0.92. 

VIII. CONCLUSIONS 

We have in detail studied the excitation spectrum for a doubly quantized vortex in a 
trapped condensate. The instability regions were studied numerically for a wide range of trap 
shapes and interaction strengths. For the previously studied two-dimensional case, which is 
expected to described pancake shaped condensates we explained the instabilities in terms 
of level crossings between the core mode and the quadrupole modes of the condensates and 
found an analytical approximation for the position of the instability regions. A corresponding 
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FIG. 11: Dynamical stability regions for a doubly quantized vortex in a two-dimensional anhar- 
monic trap, as a function of coupling strength C and anharmonicity a. 

study of the anharmonic trap was carried out in order to point out the connections and 
differences between spectral and energetic stability. It was found that the doubly quantized 
vortex is in this case stable in the weak-coupling limit, but for stronger coupling the spectrum 
is similar to that for a harmonically trapped condensate. 

We have systematically mapped out the regions of instability in a three-dimensional trap 
for a wide range of aspect ratio. In the cigar shaped regime, as the interaction strength 
becomes larger the unstable modes acquire successively more nodes in the axial direction. 
Comparison of the imaginary parts of the computed mode frequencies with the results of 
the experiment performed by Shin et al. [3], was seen to yield qualitative agreement. 
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